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MONTE  CARLO  SIMULATION  OF  LASER  BACKSCATTER 
FROM  SEA  WATER 

B.W.  Koerber  and  D.M.  Phillips 


SUMMARY 

A  Monte  Carlo  simulation  study  of  laser  backscatter  from 
sea  water  has  been  carried  out  to  provide  data  required  to 
assess  the  feasibility  of  measuring  inherent  optical 
propagation  properties  of  sea  water  from  an  aircraft.  The 
possibility  was  examined  of  deriving  such  information  from 
the  backscatter  component  of  the  return  signals  measured 
by  the  WRELADS  laser  airborne  depth  sounder  system. 
Computations  were  made  for  various  water  turbidity 
conditions  and  for  different  fields  of  view  of  the  WRELADS 
receiver.  Using  a  simple  model  fitted  to  the  computed 
backscatter  data,  it  was  shown  that  values  of  the 
scattering  and  absorption  coefficients  can  be  derived  from 
the  initial  amplitude  and  the  decay  rate  of  the 
backscatter  envelope. 
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1 .  INTRODUCTION 

Several  experimental  studies  have  recently  been  carried  out  at  the  Defence 
Research  Centre  Salisbury  to  determine  the  effects  of  water  turbidity  on  the 
performance  of  the  WRELADS  laser  airborne  depth  sounder(ref . 1-3) .  Combined 
aircraft  and  boat  trials  were  held  to  obtain  simultaneous  measurements  of 
return  signals  by  the  green  receiver  of  the  WRELADS  system  and  of  the 
attenuation  coefficient  of  the  waters  in  which  the  trials  were  conducted.  An 
empirical  relationship  was  established  between  the  maximum  depths  measured  by 
the  WRELADS  system  and  the  values  of  attenuation  coefficient  of  the  sea  water. 
Analysis  of  the  backscatter  component  of  the  green  receiver  return  signals 
indicated  that  the  amplitude  and  the  decay  rate  of  these  signals  might  be 
suitable  parameters  for  determining  water  turbidity  from  an  aircraft. 

To  provide  further  insight  into  the  feasibility  of  measuring  optical 
properties  of  sea  water  from  an  aircraft,  Monte  Carlo  simulation  studies  of 
laser  backscatter  from  sea  water  have  been  undertaken  using  a  program 
developed  at  DRCS  for  use  on  a  large  digital  computer.  Information  was 
obtained  from  this  model  of  the  magnitude  of  the  backscattered  light  signal 
versus  time  produced  at  an  airborne  receiver  from  a  pulse  of  laser  light. 
Computations  were  performed  covering  a  large  range  of  water  turbidity 
conditions  corresponding  to  the  range  in  which  WRELADS  is  operated.  Computed 
values  of  backscatter  amplitude  and  backscatter  decay  rate  were  then  compared 
with  the  values  of  the  scattering,  absorption  and  attenuation  data  used  as 
input  to  the  program.  Using  these  results,  it  was  possible  to  confirm  the 
feasibility  of  deriving  information  of  light  propagation  parameters  for  sea 
water  from  airborne  measurements  of  laser  backscatter  and  to  establish 
experimental  conditions  under  which  valid  information  can  be  obtained. 


2.  MONTE  CARLO  MODEL 
2.1  Simulation  of  WRELADS  operation 

The  computations  of  backscattered  light  from  a  laser  beam  required  for  this 
study  were  made  using  a  Monte  Carlo  model  developed  for  simulation  studies 
of  the  WRELADS  system.  This  model  has  been  programmed  for  the  IBM  370 
computer  which  provides  adequate  speed  and  storage  capacity  to  handle  the 
large  number  of  complicated  calculations  required. 

During  operations  of  the  WRELADS  equipment  (see  figure  1),  the  airborne 
laser  transmitter  and  optical  receiver  are  flown  at  a  height  of  500  m  above 
the  sea  surface.  The  narrow  transmitter  beam  consisting  of  infra  red  and 
green  components  produced  by  a  frequency  doubled  Nd:YAG  laser  is  directed 
down  at  the  sea  in  a  series  of  rapid  pulses  of  short  duration.  The  1064  nm 
infra  red  beam  is  reflected  at  the  water  surface  and  the  return  signal 
measured  by  an  infra  red  receiver,  providing  a  datum  for  the  depth  sounding 
measurement.  The  532  nm  green  beam  penetrates  the  water  medium  to  the  sea 
bottom  where  a  signal  is  reflected  back  to  the  green  receiver  which  has  an 
adjustable  field  qf  view  covering  an  area  of  up  to  10  m  radius  at  the  sea 
surface.  An  estimate  of  the  depth  of  the  sea  bottom  is  obtained  from  the 
time  interval  between  the  surface  and  bottom  reflections. 

The  return  signal  at  the  green  receiver  consists  of  four  components :  the 
surface  reflection,  backscatter,  the  sea  bottom  reflection  and  general 
background  reflected  sunlight.  It  is  the  backscatter  signal  produced  by 
multiple  scattering  in  the  water  between  the  surface  and  sea  bottom  that  is 
of  relevance  to  this  present  investigation  into  the  feasibility  of 
monitoring  optical  properties  of  water  from  an  aircraft.  The  modelling  and 
analysis  of  these  backscatter  signals  are  described  in  the  following 
sections. 
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2.2  Description  of  the  model 

For  the  purpose  of  the  modelling  (see  figure  2),  the  laser  transmitter  beam 
is  assumed  to  be  a  vertical,  parallel  light  beam  of  infinitesimal  width. 
The  laser  emission  is  modelled  as  an  instantaneous  pulse  of  light.  To 
simplify  the  analysis,  flat  sea  conditions  (no  waves  or  ripples)  are 
assumed. 

The  model  simulates  the  random  passage  in  the  water  of  individual  photons 
from  the  transmitter  light  beam  through  series  of  scatterings  until  the 
photons  reach  the  boundaries  of  the  sea  water  medium  or  move  laterally  well 
outside  the  field  of  view  of  the  airborne  optical  receiver.  For  each 
segment  of  a  photon  trajectory,  values  are  determined  for  the  distance 
travelled  before  scattering  and  for  the  scattering  angles  by  making  random 
selections  from  cumulative  probability  data  for  these  parameters.  Random 
numbers  required  for  this  procedure  are  generated  using  the  multiplicative 
congruential  method  in  a  subroutine  of  the  program.  Weighting  factors  are 
computed  for  each  photon  trajectory  to  account  for  the  absorption  losses 
along  that  path.  Transit  times  through  the  water  are  computed  for  each 
photon  path.  After  the  computations  have  been  repeated  for  a  large  number 
of  photons,  a  space-time  map  is  built  up  of  the  distribution  of  light  flux 
in  the  water. 

For  each  photon  that  is  backscattered  from  the  sea  water,  tests  are  carried 
out  in  the  program  to  determine  whether  this  photon  is  within  the  field  of 
view  and  travelling  towards  the  entrance  aperture  of  the  airborne  receiver. 
Allowance  is  made  for  a  change  in  direction  at  the  sea  surface  due  to 
refraction  at  this  interface.  Totals  are  kept  of  the  numbers  of  photons 
(weighted  for  the  effects  of  absorption)  collected  at  the  receiver  during 
successive  5  ns  intervals  of  photon  travel  time.  To  examine  the  effect  of 
changing  the  field  of  view  of  the  receiver,  sub- totals  are  kept  of  the 
numbers  of  these  photons  collected  from  concentric  annular  regions  at  the 
sea  surface  having  radii  0.25,  0.50,  1.0,  2,0,  5.0  and  10.0  m. 

Very  large  numbers  of  photon  trajectories  must  be  computed  to  reduce  the 
chance  variations  produced  in  the  results  due  to  the  random  Monte  Carlo 
process  used  in  this  simulation  study.  In  order  to  keep  computing  times 
within  manageable  limits,  it  was  decided  to  substantially  increase  the 
effective  collecting  area  of  the  optical  receiver  in  the  model,  thereby 
producing  considerable  smoothing  of  results  by  averaging  over  the  much 
larger  area.  Hence  in  the  model,  the  photons  are  counted  as  they  arrive  in 
the  concentric  annular  regions  in  the  receiver  aperture  plane  having 
radii  10,  20,  30,  40  and  50  m.  Totals  obtained  for  these  larger  areas  can 
be  scaled  to  provide  results  relevant  to  the  actual  receiver  size  (180  mm 
diameter) . 

The  computer  program  is  executed  in  a  series  of  separate  runs  and  after 
each  run  the  output  data  are  combined  in  a  single  data  set.  As  each 
additional  run  is  carried  out,  an  assessment  is  made  of  the  magnitude  of 
the  random  errors  in  the  data  due  to  the  Monte  Carlo  process.  These  errors 
are  progressively  reduced  as  the  data  sample  size  is  increased.  The 
computations  are  continued  until  the  random  errors  are  reduced  to  an 
acceptable  level. 


3.  CALCULATIONS  OF  LASER  BACKSCATTER 
3.1  Data  on  optical  properties  of  water 

The  computations  of  backscattered  light  from  a  laser  beam  were  carried  out 
for  i  range  of  water  turbidity  conditions  from  clear  water  to  the  highest 
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turbidity  in  which  WRELADS  can  be  operated  effectively.  Computations  were 
performed  for  values  of  attenuation  coefficient,  c,  of  0.1,  0.5,  2.0  and 
5.0  m-1. 

The  magnitude  of  the  contributions  of  absorption  and  scattering  to  these 
levels  of  attenuation  were  determined  by  referring  to  an  experimental  study 
of  absorption  versus  attenuation  carried  out  by  DRCS  during  1979  and  1980. 
In  this  study,  a  series  of  in  situ  measurements  were  made  at  several 
wavelengths  of  the  absorption  coefficient,  a,  and  the  attenuation 
coefficient,  c,  in  waters  of  various  turbidity  in  the  Gulf  of  St  Vincent. 
The  measurements  were  made  from  a  small  boat  using  an  absorption  meter 
instrument  (figure  3)  and  a  transmissometer  instrument  (figure  4,  and  also 
see  reference  4).  Both  instruments  were  developed  and  made  at  DRCS:  the 
absorption  meter  is  based  on  a  design  produced  at  the  Stanford  Research 
Institute(ref .5)  and  the  transmissometer  on  a  design  by  Austin  and  Petzold 
at  the  Scripps  Institution  of  0ceanography(ref .6) .  Data  of  a  versus  c 
measured  at  a  wavelength  of  530  nm  are  shown  in  figure  5.  For  comparison, 
several  examples  of  data  measured  by  the  Scripps  Institution  of 
Oceanography (ref .7)  are  included  in  the  figure;  there  is  general  agreement 
between  the  two  sets  of  data,  but  the  values  of  a  appear  to  be  somewhat 
higher  at  small  c  in  the  Scripps'  data.  The  values  of  a  used  in  the  Monte 
Carlo  study  were  obtained  from  the  regression  line  fitted  to  the  DRCS  data. 
Values  of  the  scattering  coefficient,  b,  were  then  calculated  by 
subtracting  the  values  of  a  from  the  corresponding  values  of  c.  Table  1 
provides  a  list  of  the  values  of  c,  a  and  b  used  in  the  Monte  Carlo  study. 


TABLE  1.  INPUT  DATA  ON  OPTICAL  PROPERTIES  OF  WATER 


Attenuation 
Coefficient 
c  (m-1) 

Absorption 
Coefficient 
a  (m-1) 

Scattering 
Coefficient 
b  (m-1) 

0.1 

0.0725 

0.0275 

0.5 

0.128 

0.372 

2.0 

0.337 

1.663 

5.0 

0.754 

4.246 

Information  of  the  volume  scattering  function,  0(0),  was  obtained  from 
experimental  data  published  by  the  Scripps  Institution  of 
Oceanography (ref . 7) .  In  the  model,  basic  $(6)  data  were  provided  in  the 
form  of  tables  for  seven  representative  curves  measured  at  a  wavelength  of 
530  nm  in  waters  having  a  range  of  values  of  b  from  0.009  to  1.818  m-1. 
During  the  operation  of  the  model,  appropriate  data  of  0(0)  were  generated 
for  the  particular  values  of  b  required  in  the  simulation  study  by 
interpolation  (or  extrapolation)  of  the  tabulated  basic  data.  The  accuracy 
of  this  operation  was  significantly  improved  by  first  normalizing  the  0(0) 
data  by  dividing  by  the  corresponding  values  of  b  to  produce  nearly 
constant-valued  functions  versus  b  for  each  of  the  tabulated  values  of  the 
scattering  angle,  8  (examples  are  shown  in  figure  6).  Values  of  cumulative 
probability,  P(0),  were  then  obtained  from  the  derived  0(0)  data  using  the 
trapezoid  integration  technique  to  evaluate  the  following  equation: 


2n/Jp(6)  sin0  d0 
2n/Jp(0)  sin©  d0 


2n/Jp(0)  sin©  d0 

b  (1) 


The  values  of  P(0)  obtained  for  b  =  1.663  m-1  are  shown  in  figure  7. 

3.2  Analysis  of  computed  backscatter  data 

The  amplitude  of  the  backscattered  light  signal  decays  exponentially  with 
time  measured  from  the  entry  of  the  laser  beam  into  the  water.  In  clear 
water,  the  initial  amplitude  of  the  signal  is  small  and  the  decay  time  is 
long,  but  in  turbid  water,  the  initial  amplitude  is  large  and  the  decay 
time  is  short.  In  a  previous  study  of  the  backscatter  signals  measured  by 
the  WRELADS  system(ref .2) ,  a  simple  model  had  been  suggested  to  describe 
the  shape  of  the  backscatter  signal.  The  amplitude  of  the  backscattered 
light  is  given  by: 


Bz  =  B0  exp  (-21^)  (2) 

where  BQ  is  the  initial  backscatter  signal, 

k k  is  the  backscatter  attenuation  coefficient, 

and  Z  is  the  depth  below  the  water  surface. 

Signals  recorded  as  a  function  of  time  can  be  tranformed  to  a  function  of 
depth  using  the  value  of  the  velocity  of  light  in  water  (2.245  x  108  ms-1) 
as  the  conversion  factor. 

The  purpose  of  the  Monte  Carlo  simulation  study  of  laser  backscatter  was  to 
provide  data  against  which  the  validity  of  the  above  simple  model  (or  a 
modified  form  of  it)  could  be  assessed.  It  was  expected  that  the 
simulation  study  results  would  provide  a  means  of  establishing 
relationships  between  the  model  parameters  B^  and  k^  and  the  parameters 

describing  the  inherent  optical  properties  of  the  water  viz,  a,  b,  c,  and 
0(0).  As  a  step  towards  achieving  this  aim,  the  simple  backscatter  model 
was  developed  a  stage  further.  Assuming  the  backscattered  light 
originating  from  just  below  the  water  surface  is  essentially  light  that  has 
been  scattered  once  only,  theory  indicates  that  the  parameter  B^  would  be 

proportional  to  0(180°)  and  hence  0(180°)  can  be  evaluated  from  Bq.  Values 

of  0(180°)  can  in  turn  be  converted  to  values  of  b  using  an  empirical  (non¬ 
linear)  relationship  between  these  parameters  derived  from  the  published 
Scripps'  data  (figure  8).  Thus  the  model  was  extended  to  provide 
information  of  b  from  the  values  of  obtained  by  fitting  the  model  to  the 

backscatter  data.  Any  significant  deviations  of  these  estimated  values 
of  b  from  the  actual  values  used  as  input  to  the  Monte  Carlo  model  could 
then  be  attributed  to  a  breakdown  of  the  simple  model  due  to  the  effects  of 
multiple  scattering.  Fitting  the  model  to  the  Monte  Carlo  backscatter  data 
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was  carried  out  by  first  taking  logarithms  of  equation  2: 


in  Bz  =  in  Bq  -  21^ 


(3) 


Linear  regression  analysis  was  then  used  on  the  backscatter  data  computed 
for  5  ns  intervals  during  the  initial  100  ns  of  photon  travel  time 
(corresponding  to  a  depth  of  11.23  m) .  For  each  time  interval,  an 
'effective  depth'  was  calculated  allowing  for  the  exponential  decrease  of 
light  flux  collected  during  the  interval.  From  the  fitted  regression  line, 
values  were  obtained  of  k^  and  In  B^  and  hence  of  b.  The  standard  errors  in 
these  parameters  were  computed  Iron)  the  magnitude  of  the  spread  of  the  data 
points  about  the  regression  line  (this  spread  is  due  mainly  to  the  random 
Monte  Carlo  process).  The  analysis  procedure  was  carried  out  on  the  data 
obtained  for  each  of  the  six  different  receiver  fields  of  view  for  which 
the  Monte  Carlo  simulations  had  been  performed  (for  the  two  smallest  field 
of  view  settings,  only  data  collected  in  the  innermost  receiver  annuli  were 
used  to  avoid  'vignetting  effects'  on  the  backscattered  light  produced  at 
the  greater  depths. 


4.  RESULTS 

The  Monte  Carlo  computations  were  continued  until  the  random  errors  in  the 
results  were  reduced  to  an  acceptable  level.  It  was  found  that  these  errors 
had  levelled  off  after  12  runs,  each  of  duration  of  approximately  10  min  CPU 
time.  The  final  errors  in  the  estimated  values  of  a  and  b  obtained  for  the 
maximum  receiver  field  of  view  (10  m  radius)  were  about  3  and  12%  respectively 
for  c  =  0.5,  2.0  and  5.0  m-1.  The  errors  were  approximately  doubled  for 
c  =  0.1  m-1.  The  total  numbers  of  photon  trajectories  computed  were  16.8,  3.0, 
0.6  and  0.6  millions  respectively,  for  c  =  0.1,  0.5,  2.0  and  5.0  m-1. 

Figure  9  shows  the  numbers  of  backscattered  photons  versus  time  collected  at 
the  airborne  receiver.  Data  are  plotted  for  the  maximum  receiver  field  of  view 
(10  m  radius).  Also  shown  on  the  figures  are  the  regression  lines  fitted  to 
the  backscatter  data. 

Computed  values  of  the  backscatter  attenuation  coefficient,  k^,  versus  radius 

of  the  receiver  field  of  view  are  shown  in  figures  10  to  13.  The  values  of  k^ 

are  seen  to  lie  between  the  values  of  the  absorption  coefficient,  a,  and  the 
attenuation  coefficient,  c.  As  the  field  of  view  of  the  receiver  is  increased, 
the  value  of  k^  progressively  decreases  and  tends  to  the  value  of  the 

absorption  coefficient  at  the  maximum  field  of  view  setting  (10  m  radius). 
The  trend  shown  in  these  data  is  statistically  significant,  as  is  indicated  by 
the  relatively  small  standard  errors  (shown  as  error  bars  in  the  figures) 
applying  to  the  data. 

The  computed  values  of  the  scattering  coefficient,  b,  versus  radius  of  the 
receiver  field  of  view  are  shown  in  figure  14.  These  values  are  distributed  on 
both  sides  of  the  actual  values  of  b  used  as  input  to  the  Monte  Carlo  model. 
The  standard  errors  in  the  computed  values  are  of  the  same  magnitude  as  the 
differences  between  the  computed  and  actual  values.  Therefore,  it  is  inferred 
that  the  apparent  variations  are  due  to  chance  and  the  computed  values  are 
valid  estimates  of  the  scattering  coefficient  (the  source  of  these  chance 
variations  is  largely  from  the  random  Monte  Carlo  process).  However,  this 
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general  conclusion  possibly  may  not  be  valid  for  the  data  obtained  for 
c  =  5.0  m-1.  The  computed  scattering  coefficient  data  for  this  turbidity  level 
appear  to  have  a  significant  bias  towards  higher  values  than  the  actual  input 
value.  This  may  be  evidence  that  the  simple  model  of  backscatter  does 
breakdown  at  very  high  water  turbidity  levels  due  to  the  effects  of  multiple 
scattering.  Nevertheless,  the  model  appears  to  be  adequate  over  virtually  the 
entire  range  of  water  turbidity  conditions  in  which  the  WRELADS  system  is 
normally  operated. 


5.  CONCLUSIONS 

The  following  conclusions  are  obtained  from  the  results  of  the  Monte  Carlo 
simulation  study  of  backscatter  from  a  laser  beam. 

(1)  It  is  feasible  to  monitor  optical  propagation  properties  of  sea  water 
by  making  measurements  of  laser  backscatter  from  an  aircraft.  The  WRELADS 
system  is  suitable  for  such  measurements. 

(2)  The  simple  model  of  backscatter  examined  in  this  study  is  adequate  to 
describe  the  shape  of  the  backscatter  signal  over  the  range  of  water 
turbidity  conditions  in  which  WRELADS  is  operated. 

(3)  Accurate  estimates  of  the  absorption  coefficient,  a,  can  be  obtained 
from  measurements  of  the  decay  rate  of  the  backscatter  envelope  providing 
an  adequately  large  receiver  field  of  view  is  used.  For  the  WRELADS 
receiver  operating  at  a  height  of  500  m,  the  greatest  accuracy  is  obtained 
using  the  maximum  field  of  view  of  radius  10  m  at  the  sea  surface. 

(4)  Measurements  of  the  initial  magnitude  of  the  backscatter  signal  can  be 
used  to  obtain  estimates  of  the  scattering  properties  of  sea  water. 
Accurate  estimates  can  be  directly  obtained  of  the  volume  scattering 
function  at  0  =  180°.  These  in  turn  can  be  converted  to  values  of  the 
scattering  coefficient,  b,  using  an  empirical  (non-linear)  relationship 
between  these  parameters.  The  accuracy  of  the  estimates  of  b  will  depend, 
of  course,  on  the  general  applicability  of  the  empirical  relationship  used 
for  the  conversions.  This  is  an  area  where  more  experimental  data  would  be 
beneficial . 

(5)  Further  experimental  studies  should  be  carried  out  using  the  WRELADS 
system  to  develop  the  measurement  techniques  required  for  routine 
monitoring  of  optical  properties  of  sea  water  from  an  aircraft.  An 
assessment  is  required  of  the  accuracy  of  data  obtained  under  normal 
operating  conditions. 
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Figure  2 


Figure  2.  Geometry  used  for  Monte  Carlo  simulation  model 
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CUMULATIVE  PROBABILITY,  P(0) 


SCATTERING  COEFFICIENT,  b  (m-') 


Figure  8.  Relationship  between  values  of  volume  scattering 
function  at  180°  and  scattering  coefficient 
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Figure  10.  Values  of  backscatter  attenuation 

coefficient  computed  for  c  =  0.1  m  1 
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Figure  14.  Values  of  scattering  coefficient  computed  for 
c  =  0.1,  0.5,  2.0  and  5.0m  1 
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